Getting ready
library(dada2); packageVersion("dada2")
[1] ‘1.0.3’
library(phyloseq); packageVersion("phyloseq")
[1] ‘1.16.2’
library(stringr); packageVersion("stringr")
[1] ‘1.0.0’
library(ggplot2); packageVersion("ggplot2")
[1] ‘2.1.0’
library(dplyr); packageVersion("dplyr")
[1] ‘0.5.0’
Load data from Nature 488, pp. 621-626: STAT.RData
path <- "../data/STAT.rdata"
load(path)
phy
phyloseq-class experiment-level object
otu_table() OTU Table: [ 6547 taxa and 96 samples ]
sample_data() Sample Data: [ 96 samples by 6 sample variables ]
tax_table() Taxonomy Table: [ 6547 taxa by 9 taxonomic ranks ]
phy_tree() Phylogenetic Tree: [ 6547 tips and 6321 internal nodes ]
head(phenotypes)
GIP Insulin Leptin IGF1 BMD mFat pFat
C1 8.61 1 2320 1250 0.0492 4.0 19.7
C2 28.60 137 1040 2180 0.0487 3.7 18.6
C3 19.20 543 1670 2160 0.0510 3.6 19.5
C4 38.30 606 972 2070 0.0465 4.0 21.3
C5 13.20 1 2210 1240 0.0547 4.1 21.6
C6 36.80 696 2330 1830 0.0512 3.7 17.7
tail(phenotypes)
GIP Insulin Leptin IGF1 BMD mFat pFat
VP5 49.5 418 2130 1810 0.0484 4.2 21.9
VP6 49.9 271 1740 1610 0.0476 4.7 23.9
VP7 48.6 224 1710 1520 0.0482 5.5 26.1
VP8 37.5 1 795 2020 0.0504 4.0 19.9
VP9 41.7 1200 4310 1390 0.0484 4.8 27.6
VP10 49.1 1 917 1970 0.0509 5.3 25.5
head(sample_data(phy))
Sample Data: [6 samples by 6 sample variables]:
X.SampleID BarcodeSequence LinkerPrimerSequence Treatment Location
cecal_C1 cecal_C1 NA NA C cecal
cecal_C10 cecal_C10 NA NA C cecal
cecal_C2 cecal_C2 NA NA C cecal
cecal_C3 cecal_C3 NA NA C cecal
cecal_C4 cecal_C4 NA NA C cecal
cecal_C5 cecal_C5 NA NA C cecal
Description
cecal_C1 missing_description
cecal_C10 missing_description
cecal_C2 missing_description
cecal_C3 missing_description
cecal_C4 missing_description
cecal_C5 missing_description
levels(sample_data(phy)$Treatment)
[1] "C" "P" "T" "V" "VP"
Merge phenotypes variable with the phyloseq data
unique(sample_names(phy))
[1] "cecal_C1" "cecal_C10" "cecal_C2" "cecal_C3" "cecal_C4" "cecal_C5"
[7] "cecal_C6" "cecal_C7" "cecal_C8" "cecal_C9" "cecal_P1" "cecal_P10"
[13] "cecal_P2" "cecal_P3" "cecal_P4" "cecal_P5" "cecal_P6" "cecal_P7"
[19] "cecal_P8" "cecal_P9" "cecal_T1" "cecal_T10" "cecal_T2" "cecal_T3"
[25] "cecal_T4" "cecal_T5" "cecal_T6" "cecal_T7" "cecal_T8" "cecal_T9"
[31] "cecal_V1" "cecal_V10" "cecal_V2" "cecal_V3" "cecal_V4" "cecal_V5"
[37] "cecal_V6" "cecal_V7" "cecal_V8" "cecal_V9" "cecal_VP1" "cecal_VP10"
[43] "cecal_VP2" "cecal_VP3" "cecal_VP4" "cecal_VP5" "cecal_VP6" "cecal_VP7"
[49] "cecal_VP8" "cecal_VP9" "fecal_C1" "fecal_C10" "fecal_C2" "fecal_C3"
[55] "fecal_C4" "fecal_C5" "fecal_C6" "fecal_C7" "fecal_C8" "fecal_C9"
[61] "fecal_P1" "fecal_P10" "fecal_P2" "fecal_P3" "fecal_P4" "fecal_P6"
[67] "fecal_P7" "fecal_P8" "fecal_P9" "fecal_T1" "fecal_T10" "fecal_T2"
[73] "fecal_T3" "fecal_T4" "fecal_T6" "fecal_T7" "fecal_T8" "fecal_T9"
[79] "fecal_V1" "fecal_V10" "fecal_V2" "fecal_V3" "fecal_V4" "fecal_V5"
[85] "fecal_V6" "fecal_V7" "fecal_V8" "fecal_V9" "fecal_VP1" "fecal_VP10"
[91] "fecal_VP2" "fecal_VP3" "fecal_VP4" "fecal_VP7" "fecal_VP8" "fecal_VP9"
pheno1 <- phenotypes
pheno2 <- phenotypes
rownames(pheno1) <- paste('cecal', rownames(pheno1), sep = "_")
rownames(pheno2) <- paste('fecal', rownames(pheno2), sep = "_")
pheno <- rbind(pheno1, pheno2)
head(pheno)
GIP Insulin Leptin IGF1 BMD mFat pFat
cecal_C1 8.61 1 2320 1250 0.0492 4.0 19.7
cecal_C2 28.60 137 1040 2180 0.0487 3.7 18.6
cecal_C3 19.20 543 1670 2160 0.0510 3.6 19.5
cecal_C4 38.30 606 972 2070 0.0465 4.0 21.3
cecal_C5 13.20 1 2210 1240 0.0547 4.1 21.6
cecal_C6 36.80 696 2330 1830 0.0512 3.7 17.7
tail(pheno)
GIP Insulin Leptin IGF1 BMD mFat pFat
fecal_VP5 49.5 418 2130 1810 0.0484 4.2 21.9
fecal_VP6 49.9 271 1740 1610 0.0476 4.7 23.9
fecal_VP7 48.6 224 1710 1520 0.0482 5.5 26.1
fecal_VP8 37.5 1 795 2020 0.0504 4.0 19.9
fecal_VP9 41.7 1200 4310 1390 0.0484 4.8 27.6
fecal_VP10 49.1 1 917 1970 0.0509 5.3 25.5
# phenotypes_fixed <- phenotypes
# rownames(phenotypes_fixed) <- str_replace_all(rownames(phenotypes),
# "C", "cecal_C")
# rownames(phenotypes_fixed) <- str_replace_all(rownames(phenotypes_fixed),
# "C", "cecal_C")
phy2 <- merge_phyloseq(phy, sample_data(pheno))
phy
phyloseq-class experiment-level object
otu_table() OTU Table: [ 6547 taxa and 96 samples ]
sample_data() Sample Data: [ 96 samples by 6 sample variables ]
tax_table() Taxonomy Table: [ 6547 taxa by 9 taxonomic ranks ]
phy_tree() Phylogenetic Tree: [ 6547 tips and 6321 internal nodes ]
phy2
phyloseq-class experiment-level object
otu_table() OTU Table: [ 6547 taxa and 96 samples ]
sample_data() Sample Data: [ 96 samples by 13 sample variables ]
tax_table() Taxonomy Table: [ 6547 taxa by 9 taxonomic ranks ]
phy_tree() Phylogenetic Tree: [ 6547 tips and 6321 internal nodes ]
head(sample_data(phy2))
Sample Data: [6 samples by 13 sample variables]:
X.SampleID BarcodeSequence LinkerPrimerSequence Treatment Location
cecal_C1 cecal_C1 NA NA C cecal
cecal_C10 cecal_C10 NA NA C cecal
cecal_C2 cecal_C2 NA NA C cecal
cecal_C3 cecal_C3 NA NA C cecal
cecal_C4 cecal_C4 NA NA C cecal
cecal_C5 cecal_C5 NA NA C cecal
Description GIP Insulin Leptin IGF1 BMD mFat pFat
cecal_C1 missing_description 8.61 1 2320 1250 0.0492 4.0 19.7
cecal_C10 missing_description 27.70 516 2680 3720 0.0479 4.0 19.8
cecal_C2 missing_description 28.60 137 1040 2180 0.0487 3.7 18.6
cecal_C3 missing_description 19.20 543 1670 2160 0.0510 3.6 19.5
cecal_C4 missing_description 38.30 606 972 2070 0.0465 4.0 21.3
cecal_C5 missing_description 13.20 1 2210 1240 0.0547 4.1 21.6
Estimate observed species with Chao1 and Shannon diversity using estimate_richness() function
alpha_meas <- c("Observed", "Chao1", "Shannon")
alpha_div <- estimate_richness(phy2, measures = alpha_meas)
glimpse(alpha_div)
Observations: 96
Variables: 4
$ Observed <dbl> 662, 592, 617, 731, 601, 698, 579, 529, 462, 617, 484, 593, ...
$ Chao1 <dbl> 1650.4286, 1508.4531, 1216.3939, 1746.2188, 1452.0685, 1434....
$ se.chao1 <dbl> 148.60209, 153.81056, 91.48394, 144.87205, 137.36109, 104.96...
$ Shannon <dbl> 4.836219, 4.814829, 4.775428, 4.918447, 4.636482, 4.944320, ...
Subsetting
phy3 <- subset_samples(phy2, sample_sums(phy2) > 3500)
phy3
phyloseq-class experiment-level object
otu_table() OTU Table: [ 6547 taxa and 93 samples ]
sample_data() Sample Data: [ 93 samples by 13 sample variables ]
tax_table() Taxonomy Table: [ 6547 taxa by 9 taxonomic ranks ]
phy_tree() Phylogenetic Tree: [ 6547 tips and 6321 internal nodes ]
alpha_div_sub <- estimate_richness(phy3, measures = alpha_meas)
Estimate rarified diversities in #3 averaged over 20 replicates and compare the two estimates
set.seed(0)
alpha_div_rare <- matrix(0, ncol = ncol(alpha_div), nrow = nrow(alpha_div))
for (i in 1:20) {
alpha_div_rare <- alpha_div_rare +
estimate_richness(
rarefy_even_depth(phy3, sample.size = 3500,
verbose = F, trimOTUs = F),
measures = alpha_meas
)
}
alpha_div_rare <- alpha_div_rare / 20
glimpse(alpha_div_rare)
Observations: 93
Variables: 4
$ Observed <dbl> 467.70, 409.80, 443.30, 501.20, 398.80, 481.90, 376.10, 354....
$ Chao1 <dbl> 643.5768, 550.6086, 598.5164, 705.6298, 563.2645, 668.5343, ...
$ se.chao1 <dbl> 33.85371, 29.96853, 31.54622, 37.52648, 34.22859, 34.91094, ...
$ Shannon <dbl> 4.732475, 4.723733, 4.675471, 4.806043, 4.540449, 4.826177, ...
plot(alpha_div$Observed, alpha_div_rare$Observed)
Error in xy.coords(x, y, xlabel, ylabel, log) :
'x' and 'y' lengths differ
Summarize the data at the order level
order_phy <- tax_glom(phy3, taxrank = "Order")
glimpse(order_phy)
Formal class 'phyloseq' [package "phyloseq"] with 5 slots
..@ otu_table:Formal class 'otu_table' [package "phyloseq"] with 2 slots
..@ tax_table:Formal class 'taxonomyTable' [package "phyloseq"] with 1 slot
..@ sam_data :'data.frame': 93 obs. of 13 variables:
Formal class 'sample_data' [package "phyloseq"] with 4 slots
..@ phy_tree :List of 5
.. ..$ edge : int [1:39, 1:2] 22 22 23 24 24 25 26 26 27 27 ...
.. ..$ Nnode : int 19
.. ..$ tip.label : chr [1:21] "1948" "1931" "4911" "4220" ...
.. ..$ edge.length: num [1:39] 0.2476 0.2294 0.1725 0.4988 0.0353 ...
.. ..$ node.label : chr [1:19] "" "0.783" "0.768" "0.334" ...
.. ..- attr(*, "class")= chr "phylo"
.. ..- attr(*, "order")= chr "cladewise"
..@ refseq : NULL
Normalize the data using CLR, relative abundance, and DESeq2
order_rel <- transform_sample_counts(order_phy,
function(x) x / sum(x))
order_clr <- transform_sample_counts(order_phy,
function(x) {
y = log(x + 1)
y / sum(y)
})
glimpse(order_deseq)
Formal class 'phyloseq' [package "phyloseq"] with 5 slots
..@ otu_table:Formal class 'otu_table' [package "phyloseq"] with 2 slots
..@ tax_table:Formal class 'taxonomyTable' [package "phyloseq"] with 1 slot
..@ sam_data :'data.frame': 93 obs. of 13 variables:
Formal class 'sample_data' [package "phyloseq"] with 4 slots
..@ phy_tree :List of 5
.. ..$ edge : int [1:39, 1:2] 22 22 23 24 24 25 26 26 27 27 ...
.. ..$ Nnode : int 19
.. ..$ tip.label : chr [1:21] "1948" "1931" "4911" "4220" ...
.. ..$ edge.length: num [1:39] 0.2476 0.2294 0.1725 0.4988 0.0353 ...
.. ..$ node.label : chr [1:19] "" "0.783" "0.768" "0.334" ...
.. ..- attr(*, "class")= chr "phylo"
.. ..- attr(*, "order")= chr "cladewise"
..@ refseq : NULL
Compute appropriate univariate tests on normalized data
Location = sample_data(order_rel)$Location
rel.p = apply(otu_table(order_rel),1, function(x) wilcox.test(c(x)~Location)$p.value)
clr.p = apply(otu_table(order_clr),1, function(x) t.test(c(x)~Location)$p.value)
deseq.p = apply(otu_table(order_deseq),1, function(x) t.test(log10(c(x))~Location)$p.value)
univ.location.res = data.frame(taxa = tax_table(order_clr)[,"Order"], rel.p, clr.p, deseq.p)
head(univ.location.res)
Adjust the p-values using False Discovery Rate
univ.location.res$rel.fdr = p.adjust(univ.location.res$rel.p, method="fdr")
univ.location.res$clr.fdr = p.adjust(univ.location.res$clr.p, method="fdr")
univ.location.res$deseq.fdr = p.adjust(univ.location.res$deseq.p, method="fdr")
colSums(univ.location.res<0.05)
Warning in Ops.factor(left, right) : ‘<’ not meaningful for factors
Order rel.p clr.p deseq.p rel.fdr clr.fdr deseq.fdr
NA 10 8 9 9 7 7
univ.location.res
Order rel.p clr.p deseq.p rel.fdr
1948 Bacteroidales 1.453797e-11 9.983697e-01 4.902762e-01 8.177357e-11
1931 "Erysipelotrichales" 8.437764e-01 2.076959e-01 5.187488e-01 8.859653e-01
4911 Anaeroplasmatales 9.721986e-10 2.249341e-15 5.735335e-14 4.083234e-09
4220 Clostridiales 1.186926e-19 1.894648e-32 3.651509e-18 2.492544e-18
1026 Bacillales 8.032710e-08 1.301173e-06 4.190743e-06 2.811448e-07
1673 "Lactobacillales" 6.834879e-16 4.321846e-19 1.071339e-16 7.176623e-15
4789 Actinobacteridae 1.132115e-01 7.601720e-02 4.509191e-02 1.828801e-01
5995 Burkholderiales 7.245135e-02 8.371138e-02 1.021904e-01 1.269279e-01
5975 Neisseriales 3.437768e-01 3.224313e-01 1.034844e-01 4.010729e-01
5242 Sphingomonadales 7.253024e-02 9.877532e-02 2.072876e-01 1.269279e-01
1453 Rhodobacterales 3.437768e-01 3.224313e-01 1.066468e-01 4.010729e-01
455 Rhizobiales 3.118829e-01 3.227785e-01 8.148514e-01 4.010729e-01
6648 Myxococcales 3.118829e-01 3.227785e-01 8.165448e-01 4.010729e-01
3405 Enterobacteriales 2.819788e-02 3.464769e-02 3.530636e-02 5.921555e-02
2012 Pasteurellales 9.633972e-01 9.332040e-01 6.386468e-01 9.633972e-01
1847 Pseudomonadales 1.239806e-03 4.114033e-03 4.024734e-03 3.254492e-03
2172 Chloroplast 8.788168e-05 3.013895e-04 2.380779e-04 2.636450e-04
3806 Flavobacteriales 1.461646e-01 1.886147e-01 3.350350e-01 2.192469e-01
4609 Mycoplasmatales 5.427940e-01 2.501740e-01 2.868498e-01 5.999302e-01
3348 Desulfovibrionales 1.954624e-02 9.799033e-01 2.337175e-01 4.560790e-02
2626 Coriobacteridae 1.557592e-11 1.459745e-12 3.481561e-12 8.177357e-11
clr.fdr deseq.fdr
1948 9.983697e-01 6.052069e-01
1931 3.355087e-01 6.052069e-01
4911 1.574539e-14 4.014735e-13
4220 3.978760e-31 7.668168e-17
1026 5.464927e-06 1.760112e-05
1673 4.537938e-18 1.124905e-15
4789 1.757939e-01 1.052145e-01
5995 1.757939e-01 1.866319e-01
5975 3.765749e-01 1.866319e-01
5242 1.885711e-01 3.348492e-01
1453 3.765749e-01 1.866319e-01
455 3.765749e-01 8.165448e-01
6648 3.765749e-01 8.165448e-01
3405 9.095018e-02 9.267919e-02
2012 9.983697e-01 7.058727e-01
1847 1.234210e-02 1.207420e-02
2172 1.054863e-03 8.332725e-04
3806 3.300758e-01 4.397334e-01
4609 3.752611e-01 4.015897e-01
3348 9.983697e-01 3.505762e-01
2626 7.663660e-12 1.827820e-11
Plot the abundances using plot_heatmap(), plot_bar(), or produce box and whisker plots
plot_heatmap(order_rel, taxa.label = "Order", sample.order = "Location")+
facet_grid(.~Location, scales = "free_x")

plot_heatmap(order_clr, taxa.label = "Order", sample.order = "Location")+
facet_grid(.~Location, scales = "free_x")

plot_heatmap(order_deseq, taxa.label = "Order", sample.order = "Location")+
facet_grid(.~Location, scales = "free_x")

plot_bar(order_rel, fill="Location") +
facet_wrap(facets = "Order", ncol=5, nrow=5, scales = "free_y")

plot_bar(order_clr, fill="Location") +
facet_wrap(facets = "Order", ncol=5, nrow=5, scales = "free_y")

plot_bar(order_deseq, fill="Location") +
facet_wrap(facets = "Order", ncol=5, nrow=5, scales = "free_y")

ggplot(psmelt(order_rel), aes(x=Location, y=Abundance)) +
geom_boxplot(aes(fill=Location)) + geom_jitter() +
facet_wrap(facets="Order", ncol=5, nrow=5, scales="free_y")

ggplot(psmelt(order_clr), aes(x=Location, y=Abundance)) +
geom_boxplot(aes(fill=Location)) + geom_jitter() +
facet_wrap(facets="Order", ncol=5, nrow=5, scales="free_y")

ggplot(psmelt(order_deseq), aes(x=Location, y=Abundance)) +
geom_boxplot(aes(fill=Location)) + geom_jitter() +
facet_wrap(facets="Order", ncol=5, nrow=5, scales="free_y")

Subset the order level data to only fecal samples
fecal.rel = subset_samples(order_rel, Location = "fecal")
fecal.rel.subs = subset_taxa(fecal.rel, rowMeans(otu_table(fecal.rel))> 0.001)
Filter out taxa with low abundance (mean abundance < 0.1%)
fecal.rel.subs = subset_taxa(fecal.rel, rowMeans(otu_table(fecal.rel))> 0.001)
Perform univariate analysis of the fecal subset with respect to the Treatment variable
LS0tCnRpdGxlOiAiVVcgU3VtbWVyIEluc3RpdHV0ZXMgLSBJbnRyb2R1Y3Rpb24gdG8gTWV0YWdlbm9taWMgRGF0YSBBbmFseXNpcyAtIExhYiAyIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgojIEdldHRpbmcgcmVhZHkKCmBgYHtyIGxvYWRfcGFja2FnZXMsIG1lc3NhZ2U9RkFMU0V9CmxpYnJhcnkoZGFkYTIpOyBwYWNrYWdlVmVyc2lvbigiZGFkYTIiKQpsaWJyYXJ5KHBoeWxvc2VxKTsgcGFja2FnZVZlcnNpb24oInBoeWxvc2VxIikKCmxpYnJhcnkoc3RyaW5ncik7IHBhY2thZ2VWZXJzaW9uKCJzdHJpbmdyIikKbGlicmFyeShnZ3Bsb3QyKTsgcGFja2FnZVZlcnNpb24oImdncGxvdDIiKQpsaWJyYXJ5KGRwbHlyKTsgcGFja2FnZVZlcnNpb24oImRwbHlyIikKYGBgCgojIyBMb2FkIGRhdGEgZnJvbSBOYXR1cmUgNDg4LCBwcC4gNjIxLTYyNjogYFNUQVQuUkRhdGFgCgpgYGB7ciBsb2FkX2RhdGF9CnBhdGggPC0gIi4uL2RhdGEvU1RBVC5yZGF0YSIKbG9hZChwYXRoKQpwaHkKaGVhZChwaGVub3R5cGVzKQp0YWlsKHBoZW5vdHlwZXMpCmBgYAoKYGBge3IgaW5zcGVjdF9zYW1wbGVfZGF0YX0KaGVhZChzYW1wbGVfZGF0YShwaHkpKQpsZXZlbHMoc2FtcGxlX2RhdGEocGh5KSRUcmVhdG1lbnQpCmBgYAoKIyMgTWVyZ2UgYHBoZW5vdHlwZXNgIHZhcmlhYmxlIHdpdGggdGhlCWBwaHlsb3NlcWAgZGF0YQoKYGBge3IgcGhlbm99CnVuaXF1ZShzYW1wbGVfbmFtZXMocGh5KSkKCnBoZW5vMSA8LSBwaGVub3R5cGVzCnBoZW5vMiA8LSBwaGVub3R5cGVzCgpyb3duYW1lcyhwaGVubzEpIDwtIHBhc3RlKCdjZWNhbCcsIHJvd25hbWVzKHBoZW5vMSksIHNlcCA9ICJfIikKcm93bmFtZXMocGhlbm8yKSA8LSBwYXN0ZSgnZmVjYWwnLCByb3duYW1lcyhwaGVubzIpLCBzZXAgPSAiXyIpCnBoZW5vIDwtIHJiaW5kKHBoZW5vMSwgcGhlbm8yKQoKaGVhZChwaGVubykKdGFpbChwaGVubykKIyBwaGVub3R5cGVzX2ZpeGVkIDwtIHBoZW5vdHlwZXMKIyByb3duYW1lcyhwaGVub3R5cGVzX2ZpeGVkKSA8LSBzdHJfcmVwbGFjZV9hbGwocm93bmFtZXMocGhlbm90eXBlcyksIAojICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAiQyIsICJjZWNhbF9DIikKIyByb3duYW1lcyhwaGVub3R5cGVzX2ZpeGVkKSA8LSBzdHJfcmVwbGFjZV9hbGwocm93bmFtZXMocGhlbm90eXBlc19maXhlZCksIAojICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAiQyIsICJjZWNhbF9DIikKYGBgCgpgYGB7ciBwaHkyfQpwaHkyIDwtIG1lcmdlX3BoeWxvc2VxKHBoeSwgc2FtcGxlX2RhdGEocGhlbm8pKQpwaHkKcGh5MgpoZWFkKHNhbXBsZV9kYXRhKHBoeTIpKQpgYGAKCiMjIEVzdGltYXRlIG9ic2VydmVkIHNwZWNpZXMgd2l0aCBDaGFvMSBhbmQgU2hhbm5vbiBkaXZlcnNpdHkgdXNpbmcgYGVzdGltYXRlX3JpY2huZXNzKClgIGZ1bmN0aW9uCgpgYGB7ciBhbHBoYV9kaXZ9CmFscGhhX21lYXMgPC0gYygiT2JzZXJ2ZWQiLCAiQ2hhbzEiLCAiU2hhbm5vbiIpCmFscGhhX2RpdiA8LSBlc3RpbWF0ZV9yaWNobmVzcyhwaHkyLCBtZWFzdXJlcyA9IGFscGhhX21lYXMpCmdsaW1wc2UoYWxwaGFfZGl2KQpgYGAKCiMjIFN1YnNldHRpbmcKCmBgYHtyIHBoeTN9CnBoeTMgPC0gc3Vic2V0X3NhbXBsZXMocGh5Miwgc2FtcGxlX3N1bXMocGh5MikgPiAzNTAwKQpwaHkzCmBgYAoKYGBge3IgYWxwaGFfZGl2X3N1Yn0KYWxwaGFfZGl2X3N1YiA8LSBlc3RpbWF0ZV9yaWNobmVzcyhwaHkzLCBtZWFzdXJlcyA9IGFscGhhX21lYXMpCmBgYAoKCiMjIEVzdGltYXRlIHJhcmlmaWVkIGRpdmVyc2l0aWVzIGluICMzIGF2ZXJhZ2VkIG92ZXIgMjAgcmVwbGljYXRlcyBhbmQgY29tcGFyZSB0aGUgdHdvIGVzdGltYXRlcwoKYGBge3IgYWxwaGFfZGl2X3JhcmV9CnNldC5zZWVkKDApCmFscGhhX2Rpdl9yYXJlIDwtIG1hdHJpeCgwLCBuY29sID0gbmNvbChhbHBoYV9kaXYpLCBucm93ID0gbnJvdyhhbHBoYV9kaXYpKQpmb3IgKGkgaW4gMToyMCkgewogICAgYWxwaGFfZGl2X3JhcmUgPC0gYWxwaGFfZGl2X3JhcmUgKwogICAgICAgIGVzdGltYXRlX3JpY2huZXNzKAogICAgICAgICAgICByYXJlZnlfZXZlbl9kZXB0aChwaHkzLCBzYW1wbGUuc2l6ZSA9IDM1MDAsIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICB2ZXJib3NlID0gRiwgdHJpbU9UVXMgPSBGKSwKICAgICAgICAgICAgbWVhc3VyZXMgPSBhbHBoYV9tZWFzCiAgICAgICAgKQp9CmFscGhhX2Rpdl9yYXJlIDwtIGFscGhhX2Rpdl9yYXJlIC8gMjAKZ2xpbXBzZShhbHBoYV9kaXZfcmFyZSkKYGBgCgpgYGB7ciBwbG90X2NvbXBhcmV9CnBsb3QoYWxwaGFfZGl2JE9ic2VydmVkLCBhbHBoYV9kaXZfcmFyZSRPYnNlcnZlZCkKcGxvdChhbHBoYV9kaXYkQ2hhbzEsIGFscGhhX2Rpdl9yYXJlJENoYW8xKQpwbG90KGFscGhhX2RpdiRTaGFubm9uLCBhbHBoYV9kaXZfcmFyZSRTaGFubm9uKQpgYGAKCgojIyBTdW1tYXJpemUgdGhlIGRhdGEgYXQgdGhlIG9yZGVyIGxldmVsCgpgYGB7ciBvcmRlcl9waHl9Cm9yZGVyX3BoeSA8LSB0YXhfZ2xvbShwaHkzLCB0YXhyYW5rID0gIk9yZGVyIikKZ2xpbXBzZShvcmRlcl9waHkpCmBgYAoKIyMgTm9ybWFsaXplIHRoZSBkYXRhIHVzaW5nIENMUiwgcmVsYXRpdmUgYWJ1bmRhbmNlLCBhbmQgREVTZXEyCgpgYGB7ciBvcmRlcl9yZWx9Cm9yZGVyX3JlbCA8LSB0cmFuc2Zvcm1fc2FtcGxlX2NvdW50cyhvcmRlcl9waHksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBmdW5jdGlvbih4KSB4IC8gc3VtKHgpKQpgYGAKCmBgYHtyIG9yZGVyX2Nscn0Kb3JkZXJfY2xyIDwtIHRyYW5zZm9ybV9zYW1wbGVfY291bnRzKG9yZGVyX3BoeSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGZ1bmN0aW9uKHgpIHsKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB5ID0gbG9nKHggKyAxKQogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHkgLyBzdW0oeSkKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIH0pCmBgYAoKYGBge3Igb3JkZXJfZGVzZXF9CmxpYnJhcnkoREVTZXEyKQpjb3VudERhdGEgPSByb3VuZChhcyhvdHVfdGFibGUob3JkZXJfcGh5KSwgIm1hdHJpeCIpLCBkaWdpdHMgPSAwKQpjb3VudERhdGEgPSBjb3VudERhdGEgKyAxTApkZHMgPC0gREVTZXFEYXRhU2V0RnJvbU1hdHJpeChjb3VudERhdGEsIHNhbXBsZV9kYXRhKG9yZGVyX3BoeSksIGRlc2lnbiA9IH4xKQoKb3JkZXJfZGVzZXEgPSBtZXJnZV9waHlsb3NlcShvdHVfdGFibGUoY291bnRzKGVzdGltYXRlU2l6ZUZhY3RvcnMoZGRzKSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBub3JtYWxpemVkPVQpLCB0YXhhX2FyZV9yb3dzID0gVCksIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIHNhbXBsZV9kYXRhKG9yZGVyX3BoeSksIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIHBoeV90cmVlKG9yZGVyX3BoeSksIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIHRheF90YWJsZShvcmRlcl9waHkpKQpnbGltcHNlKG9yZGVyX2Rlc2VxKQpgYGAKCiMjIENvbXB1dGUgYXBwcm9wcmlhdGUgdW5pdmFyaWF0ZSB0ZXN0cyBvbiBub3JtYWxpemVkIGRhdGEKCmBgYHtyIHVuaXYubG9jYXRpb24ucmVzfQpMb2NhdGlvbiA9IHNhbXBsZV9kYXRhKG9yZGVyX3JlbCkkTG9jYXRpb24KcmVsLnAgPSBhcHBseShvdHVfdGFibGUob3JkZXJfcmVsKSwxLCBmdW5jdGlvbih4KSB3aWxjb3gudGVzdChjKHgpfkxvY2F0aW9uKSRwLnZhbHVlKQpjbHIucCA9IGFwcGx5KG90dV90YWJsZShvcmRlcl9jbHIpLDEsIGZ1bmN0aW9uKHgpIHQudGVzdChjKHgpfkxvY2F0aW9uKSRwLnZhbHVlKQpkZXNlcS5wID0gYXBwbHkob3R1X3RhYmxlKG9yZGVyX2Rlc2VxKSwxLCBmdW5jdGlvbih4KSB0LnRlc3QobG9nMTAoYyh4KSl+TG9jYXRpb24pJHAudmFsdWUpCgp1bml2LmxvY2F0aW9uLnJlcyA9IGRhdGEuZnJhbWUodGF4YSA9IHRheF90YWJsZShvcmRlcl9jbHIpWywiT3JkZXIiXSwgcmVsLnAsIGNsci5wLCBkZXNlcS5wKQpoZWFkKHVuaXYubG9jYXRpb24ucmVzKQpgYGAKCiMjIEFkanVzdCB0aGUgcC12YWx1ZXMgdXNpbmcgRmFsc2UgRGlzY292ZXJ5IFJhdGUKCmBgYHtyIGRvX2Zkcn0KdW5pdi5sb2NhdGlvbi5yZXMkcmVsLmZkciA9IHAuYWRqdXN0KHVuaXYubG9jYXRpb24ucmVzJHJlbC5wLCBtZXRob2Q9ImZkciIpCnVuaXYubG9jYXRpb24ucmVzJGNsci5mZHIgPSBwLmFkanVzdCh1bml2LmxvY2F0aW9uLnJlcyRjbHIucCwgbWV0aG9kPSJmZHIiKQp1bml2LmxvY2F0aW9uLnJlcyRkZXNlcS5mZHIgPSBwLmFkanVzdCh1bml2LmxvY2F0aW9uLnJlcyRkZXNlcS5wLCBtZXRob2Q9ImZkciIpCgpjb2xTdW1zKHVuaXYubG9jYXRpb24ucmVzPDAuMDUpCnVuaXYubG9jYXRpb24ucmVzCmBgYAoKIyMgUGxvdCB0aGUgYWJ1bmRhbmNlcyB1c2luZyBgcGxvdF9oZWF0bWFwKClgLCBgcGxvdF9iYXIoKWAsIG9yIHByb2R1Y2UgYm94IGFuZCB3aGlza2VyIHBsb3RzCgpgYGB7ciBwbG90X2hlYXRtYXBfcmVsfQpwbG90X2hlYXRtYXAob3JkZXJfcmVsLCB0YXhhLmxhYmVsID0gIk9yZGVyIiwgc2FtcGxlLm9yZGVyID0gIkxvY2F0aW9uIikrCiAgZmFjZXRfZ3JpZCgufkxvY2F0aW9uLCBzY2FsZXMgPSAiZnJlZV94IikKYGBgCgpgYGB7ciBwbG90X2hlYXRtYXBfY2xyfQpwbG90X2hlYXRtYXAob3JkZXJfY2xyLCB0YXhhLmxhYmVsID0gIk9yZGVyIiwgc2FtcGxlLm9yZGVyID0gIkxvY2F0aW9uIikrCiAgZmFjZXRfZ3JpZCgufkxvY2F0aW9uLCBzY2FsZXMgPSAiZnJlZV94IikKYGBgCgpgYGB7ciBwbG90X2hlYXRtYXBfZGVzZXF9CnBsb3RfaGVhdG1hcChvcmRlcl9kZXNlcSwgdGF4YS5sYWJlbCA9ICJPcmRlciIsIHNhbXBsZS5vcmRlciA9ICJMb2NhdGlvbiIpKwogIGZhY2V0X2dyaWQoLn5Mb2NhdGlvbiwgc2NhbGVzID0gImZyZWVfeCIpCmBgYAoKYGBge3IgcGxvdF9iYXJwbG90X3JlbCwgZmlnLmhlaWdodD02fQpwbG90X2JhcihvcmRlcl9yZWwsIGZpbGw9IkxvY2F0aW9uIikgKyAKICBmYWNldF93cmFwKGZhY2V0cyA9ICJPcmRlciIsIG5jb2w9NSwgbnJvdz01LCBzY2FsZXMgPSAiZnJlZV95IikKYGBgCgpgYGB7ciBwbG90X2JhcnBsb3RfY2xyLCBmaWcuaGVpZ2h0PTZ9CnBsb3RfYmFyKG9yZGVyX2NsciwgZmlsbD0iTG9jYXRpb24iKSArIAogIGZhY2V0X3dyYXAoZmFjZXRzID0gIk9yZGVyIiwgbmNvbD01LCBucm93PTUsIHNjYWxlcyA9ICJmcmVlX3kiKQpgYGAKCmBgYHtyIHBsb3RfYmFycGxvdF9kZXNlcSwgZmlnLmhlaWdodD02fQpwbG90X2JhcihvcmRlcl9kZXNlcSwgZmlsbD0iTG9jYXRpb24iKSArIAogIGZhY2V0X3dyYXAoZmFjZXRzID0gIk9yZGVyIiwgbmNvbD01LCBucm93PTUsIHNjYWxlcyA9ICJmcmVlX3kiKQpgYGAKCmBgYHtyIHBsb3RfYm94cGxvdF9yZWwsIGZpZy5oZWlnaHQ9Nn0KZ2dwbG90KHBzbWVsdChvcmRlcl9yZWwpLCBhZXMoeD1Mb2NhdGlvbiwgeT1BYnVuZGFuY2UpKSArIAogIGdlb21fYm94cGxvdChhZXMoZmlsbD1Mb2NhdGlvbikpICsgZ2VvbV9qaXR0ZXIoKSArIAogIGZhY2V0X3dyYXAoZmFjZXRzPSJPcmRlciIsIG5jb2w9NSwgbnJvdz01LCBzY2FsZXM9ImZyZWVfeSIpCmBgYAoKYGBge3IgcGxvdF9ib3hwbG90X2NsciwgZmlnLmhlaWdodD02fQpnZ3Bsb3QocHNtZWx0KG9yZGVyX2NsciksIGFlcyh4PUxvY2F0aW9uLCB5PUFidW5kYW5jZSkpICsgCiAgZ2VvbV9ib3hwbG90KGFlcyhmaWxsPUxvY2F0aW9uKSkgKyBnZW9tX2ppdHRlcigpICsgCiAgZmFjZXRfd3JhcChmYWNldHM9Ik9yZGVyIiwgbmNvbD01LCBucm93PTUsIHNjYWxlcz0iZnJlZV95IikKYGBgCgpgYGB7ciBwbG90X2JveHBsb3RfZGVzZXEsIGZpZy5oZWlnaHQ9Nn0KZ2dwbG90KHBzbWVsdChvcmRlcl9kZXNlcSksIGFlcyh4PUxvY2F0aW9uLCB5PUFidW5kYW5jZSkpICsgCiAgZ2VvbV9ib3hwbG90KGFlcyhmaWxsPUxvY2F0aW9uKSkgKyBnZW9tX2ppdHRlcigpICsgCiAgZmFjZXRfd3JhcChmYWNldHM9Ik9yZGVyIiwgbmNvbD01LCBucm93PTUsIHNjYWxlcz0iZnJlZV95IikKYGBgCgojIyBTdWJzZXQgdGhlIG9yZGVyIGxldmVsIGRhdGEgdG8gb25seSBmZWNhbCBzYW1wbGVzCgpgYGB7ciBmZWNhbC5yZWx9CmZlY2FsLnJlbCA9IHN1YnNldF9zYW1wbGVzKG9yZGVyX3JlbCwgTG9jYXRpb24gPSAiZmVjYWwiKQpgYGAKCiMjIEZpbHRlciBvdXQgdGF4YSB3aXRoIGxvdyBhYnVuZGFuY2UgKG1lYW4gYWJ1bmRhbmNlIDwgMC4xJSkKCmBgYHtyIGZlY2FsLmZlbC5zdWJzfQpmZWNhbC5yZWwuc3VicyA9IHN1YnNldF90YXhhKGZlY2FsLnJlbCwgcm93TWVhbnMob3R1X3RhYmxlKGZlY2FsLnJlbCkpPiAwLjAwMSkKYGBgCgojIyBQZXJmb3JtIHVuaXZhcmlhdGUgYW5hbHlzaXMgb2YgdGhlIGZlY2FsIHN1YnNldCB3aXRoIHJlc3BlY3QgdG8gdGhlIFRyZWF0bWVudCB2YXJpYWJsZQoKCgoKCg==